Dynamical BLUP modeling of reaction norm evolution, accommodating changing environments, overlapping generations, and multivariate data

Abstract For theoretical studies, reaction norm evolution in a changing environment can be modeled by means of the multivariate breeder's equation, with the reaction norm parameters treated as traits in their own right. This is, however, not a feasible approach for the use of field data, where the intercept and slope values are not available. An alternative approach is to use infinite‐dimensional characters and smooth covariance function estimates found by, e.g., random regression. This is difficult because of the need to find, for example, polynomial basis functions that fit the data reasonably well over time, and because reaction norms in multivariate cases are correlated, such that they cannot be modeled independently. Here, I present an alternative approach based on a multivariate linear mixed model of any order, with dynamical incidence and residual covariance matrices that reflect the changing environment. From such a mixed model follows a dynamical BLUP model for the estimation of the individual reaction norm parameter values at any given parent generation, and for updating of the mean reaction norm parameter values from generation to generation by means of Robertson's secondary theorem of natural selection. This will, for example, make it possible to disentangle the microevolutionary and plasticity components in climate change responses. The BLUP model incorporates the additive genetic relationship matrix in the usual way, and overlapping generations can easily be accommodated. Additive genetic and environmental model parameters are assumed to be known and constant, but it is discussed how they can be estimated by means of a prediction error method. The identifiability by the use of field or laboratory data containing environmental, phenotypic, fitness, and additive genetic relationship data is an important feature of the proposed model.


| INTRODUC TI ON
A reaction norm describes the phenotypes that a genotype can produce across a range of environments. Mean reaction norms in a population can evolve, and this evolution can be modeled by the use of several methods. In its simplest form, a mean reaction norm is characterized by an intercept value and a plasticity slope value, but it may also be natural to use models with a multiple of reaction norms, and they may be nonlinear.
A common model is the multivariate breeder's equation (Lande, 1979), where the mean reaction norm parameters may be treated as traits in their own right, as in, for example, Lande (2009).
For applications on field data from studies of wild populations, there are two problems with such models. First, the individual reaction norm parameters are not available, and second, additive genetic relationships in the population cannot be taken into account. The first problem can be solved by a linear transformation as shown in Ergon (2022a), and as used for comparison purposes in Section 2, but the second problem will still exist.
An alternative approach for reaction norm modeling is to use infinite-dimensional characters, as in a method introduced by Kirkpatrick and Heckman (1989), and as applied on reaction norms by Gomulkiewicz and Kirkpatrick (1992). It is then necessary to obtain smooth covariance function estimates (Kingsolver et al., 2001), and one method for that purpose is random regression (Shaeffer, 2004), where individual breeding values are modeled as relatively simple weighted sums of basis functions. Additive genetic relationship matrices can be included in such models (Oliviera et al., 2019). An obvious difficulty is here the need to find for example polynomial basis functions that fit the data reasonably well over time. Another difficulty is that the reaction norms in multivariate cases are correlated, such that they cannot be modeled independently.
In the present article, I introduce a multivariate modeling approach based on best linear unbiased predictions (BLUP), allowing for multiple (and potentially correlated) traits to have joint norms of reaction, which can be linear or approximated by power series (Gavrilets & Scheiner, 1993). This is a dynamical BLUP model in the sense that the incidence matrix for the random effects and the residual covariance matrix are functions of the changing environment.
I will develop the theory under the assumption that the parameters in the model are known, but as will be shown separately they can also be identified by a prediction error method, as introduced in a microevolutionary context in Ergon (2022aErgon ( , 2022b. This identification aspect is an essential prerequisite, that is, the model must be identifiable from available environmental, phenotypic, fitness, and additive genetic relationship data. BLUP based on linear mixed models with fixed and random effects are extensively used in domestic animal and plant breeding (Arnold et al., 2019;Ch. 26, Lynch & Walsh, 1998;Robinson, 1991). These methods may also be applied on wild populations (Kruuk, 2004;Nussey et al., 2007), although such uses have been criticized owing to errors in estimated variances of the random effects (Hadfield et al., 2010). An important application is the disentanglement of microevolutionary and plasticity components in for example climate change responses (Ergon, 2022a(Ergon, , 2022bMerilä & Hendry, 2014).
The basic BLUP equations were first developed in summation form (Henderson, 1950), but as done here, it is more convenient to use matrix formulations.
This article will be focused on how mean reaction norm parameter values, and thus mean phenotypic traits, evolve under the influence of environmental cues and changes in the fitness landscape. Such evolution of reaction norms and phenotypic traits seeks to maximize the mean fitness of a given population, and changes in the location of fitness peaks in the phenotypic space are therefore the driving force. I will thus study the dynamics of microevolutionary systems, mainly by the use of BLUP, but also with reference to the well-known multivariate breeder's equation (Lande, 1979). Although fitness can be defined as the long-run growth rate (Saether & Engan, 2015), or for nonoverlapping generations the expected geometric mean fitness (Autzen & Okasha, 2022), I will in simulations simply use the number of surviving descendants as a measure of individual fitness (Ch. 6, Rice, 2004).
It is well-known that BLUP underestimates the variances of the random effects in linear mixed models (Hadfield et al., 2010;Ch. 26, Lynch & Walsh, 1998). Here, I will show why and to which extent that is necessary in order to obtain the correct incremental changes in mean reaction norm parameter values from generation to generation.
I will also show that these changes may be found from Robertson's secondary theorem of natural selection (Robertson, 1966) applied on the estimated random effects. This also applies to nonplastic organisms, where the mean reaction norms degenerate into mean phenotypic trait values (Ergon, 2022c).
The dynamical BLUP model with Robertson updating of mean reaction norm parameter values, makes the use of the additive genetic relationship matrix A t in a standard way (Ch. 26, Lynch & Walsh, 1998). The theoretical treatment is limited to cases where only mean phenotypic traits are included in the fixed effects, and for simplicity, it assumes that generations are nonoverlapping. It is, however, also shown how cases with overlapping generations can be handled in a straightforward way. For cases with sexual reproduction, I assume a hypothetical single parent (mid-parent) occupying an intermediate phenotypic position between the two parents (Ch. 7, Rice, 2004).
The theoretical development will be general, that is, for any number of phenotypic traits and any number of environmental cues in the model. For clarity of presentation, however, some details will be

T A X O N O M Y C L A S S I F I C A T I O N
Quantitative genetics | 3 of 15 ERGON given for a system with only two phenotypic traits and two environmental cues. A similar limited system will also be used in simulations.
The additive genetic and phenotypic covariance matrices, G and P, are here assumed to be constant and known. They may, however, be estimated by means of a prediction error method (PEM), utilizing the information contained in environmental cues and individual phenotypic trait values over many generations, as well as fitness information (Ergon, 2022a(Ergon, , 2022b. With plastic traits in the dynamical BLUP model, restricted maximum likelihood (REML) methods applied on data from a single generation cannot be used for this purpose. The simple reason for this is that each element in the residual covariance matrix is a function of several nonadditive effects, such that the REML equations become indeterminate.
As developed theoretically, and verified in simulations, the dynamical BLUP model with an additive genetic relationship matrix equal to an identity matrix, that is, with random mating in an unbred population, will give the same results as a selection gradient prediction method (GRAD) based on the multivariate breeder's equation (Ergon, 2022a(Ergon, , 2022b. For large populations, these results will asymptotically also be the same as from the multivariate breeder's equation directly.
After this introduction, Theory and Methods follow in Section 2, Simulations in Section 3, and Summary and Discussion in Section 4.
Proofs of two theorems are given in Appendices A and B. For the interested reader, a user guide is given in Appendix C, including the procedure for PEM system identification. MATLAB code for the simulations is given in Supporting information.

| Notation
Mathematical symbols with descriptions in the order they appear in equations are shown in Table 1.

| Introductory example
For a simple toy example, intended to ease the readers into the concepts used below, consider a single trait y i,t measured on a single individual. In this case, the BLUP estimate of the true additive genetic value for that individual is â� i,t = h 2 y i,t − mean , where h 2 is the heritability, while mean denotes any fixed effects adjustment. If we substitute this estimate into Robertson's secondary theorem of natural selection (Ch. 6, Walsh & Lynch, 2018), we find the between- where w i,t is the relative fitness, while S is the Robertson-Price withingeneration change in the mean. In this way, we recover the standard univariate breeder's equation.
What is done below is to consider a much more complicated phenotype (an observed vector y t of individual focal traits with reaction norms) and use BLUP to estimate the vector of additive effects associated with the norm of reaction functions, with these BLUPs then substituted into the expression for Robertson's secondary theorem of natural selection. The incremental changes in the mean reaction norm parameter values thus follow from cov w i,t ,ẑ ′ i,t , where ẑ′ i,t stands for the BLUP estimates of the true additive genetic values involved (Equation (11) below).
For the special case with an additive genetic relationship matrix A t = I n , the BLUP estimates yield a matrix-based inheritance expression (using the correlated nature of the random effects) to replace h 2 , and a Robertson-Price term cov w i,t , y i,t to measure phenotypic selection (Equation (12) below). From a practical point of view, Equation (12) is unnecessary, but it is included for the purpose of comparisons with results from the multivariate breeder's equation.
TA B L E 1 Mathematical symbols with description.

Symbol Description
Δ z t Incremental change in mean trait from generation t to generation t + 1 Random effects in linear mixed model, for special case with 2 traits and 2 environmental variables,

| Background theory
For the development of the dynamical BLUP matrix equation that follows, we need some background theory. First, the Price equation for selection in a population with n individuals says that the evolution of the mean trait of an n × 1 vector z t of individual quantitative traits is described by (Price, 1970(Price, , 1972 where Δ z t = z t+1 − z t is the incremental change in mean trait value from generation to generation and where z i,t is an individual trait. Here, w i,t is the relative individual fitness, that is, individual fitness divided by the mean fitness in the population, while z descendants i,t is the mean trait of the descendants (and the parent if it survives) of individual i in generation t. The trait z i,t may be any property we can assign a numerical value to, not necessarily biological. In a biological context, the trait may be a behavioral, morphological, or physiological characteristic, but it may also be a parameter in a reaction norm model that describes a plastic organism. Disregarding the second term on the righthand side of Equation (1), we find the Robertson-Price identity, Δ z t = cov w i,t , z i,t (Robertson, 1966; Ch. 6, Walsh & Lynch, 2018), as referred to above, and which we will use below.
Second, we need to see how the multivariate breeder's equation (Lande, 1979;Lande & Arnold, 1983), where t is the selection gradient, which can be applied on the parameters in a reaction norm model. Equation 1. The vector z i,t of individual phenotypic traits is the sum of independent additive genetic effects x i,t and nonadditive environmental and genetic effects e i,t , that is, 2. The nonadditive effects e i,t are zero mean, independent, and identically distributed (iid) random variables.
3. There are no expected fitness-weighted changes in the individual additive genetic effects x i,t from one generation to the next be- influence individual fitness only through z i,t .
6. All individuals in the population are genetically unrelated, which means that the additive genetic relationship matrix A t is a unity matrix.
In what follows, we will make the use of Assumptions 1, 2, and 3, while Assumptions 4, 5, and 6 will be used only indirectly when the BLUP results with A t = I n are compared with results based on the multivariate breeder's equation.
In order to see how Equation (2) can be applied on the parameters in a reaction norm model, we may use an individual interceptslope model based on Assumptions 1 and 2 above (Lande, 2009) where u t is an environmental cue, and where the mean reaction norm iid and zero mean random variables, with variances 2 and 2 , respectively. We thus use a � i,t + i,t and b � i,t + i,t to denote individual deviations from mean values a t and b t , respectively, where a ′ i,t and b ′ i,t are the additive genetic components of these deviations. Such additive genetic deviations will be the random effects in the linear mixed model developed below, and thus the random effects that are estimated by the use of the BLUP equations.
When the reaction norm parameters in Equation (3) are treated as quantitative traits in their own right, Equation (2) leads to with G and P given by The additive genetic and phenotypic covariance matrices G and P may be time-varying, but for simplicity, we will here assume that they are constant. As discussed in Ergon (2022a), it is essential that the environmental input in Equation (3)  where selection with respect to a i,t and b i,t , as in Equation (4), is replaced by selection with respect to y i,t . Here, P yy = P aa + 2G ab u t + P bb u 2 t , while P −1 yy cov w i,t , y i,t is the selection gradient. It is essential to note that Equations (4) and (5) give identical results only asymptotically, when the population size n → ∞, and the reason for that is the differences in how the covariance functions are used in the two equations. As we will see, extensions of Equation (5) are possible for more complex reaction norm models. This is interesting because with an additive genetic relationship matrix A t = I n , the dynamical BLUP model we will develop results in incremental changes in mean reaction norm parameter values, that are identical to those found from an extended version of Equation (5). Here, we should finally note that since Equation (5) is derived from the multivariate breeder's equation (2), it is valid only under Assumptions 1-6 above.

| Development of the dynamical BLUP model
For clarity of presentation, some details will here be limited to a system with p = 2 phenotypic traits, and q = 2 environmental cues, and a similar simplified system will also be used in the simulations. The theory will, however, be developed in such a way that extensions to higher dimensions are obvious.
For a specific trait j, the individual reaction norm model with With p = 2 traits and q = 2 environmental cues, we thus obtain the linear mixed model in general form (Ch. 26, Lynch & Walsh, 1998; with y t as fixed effects vector) where e 1,t = v 1,t + 11,t u 1,t + 12,t u 2,t and , and where ⨂ is the Kronecker product operator, which means that all elements in Z t should be multiplied by I n . For use in Appendix B, we may note that It is an essential feature of the model in Equation (7) Lynch & Walsh, 1998). We thus also have that E y j,i,t = y j,t .
The random effects in Equation (7) may all be correlated, with an The residuals e 1,t and e 2,t in Equation (7) are assumed to be uncorrelated, and 2 j2 as the variances of the nonadditive effects v j,i,t , j1,i,t and j2,i,t according to Equation (6a). Also G may in general be a function of time, but we will here assume that it is constant. We will assume that G and R t are known, although they may in practice be estimated by the use of restricted maximum likelihood (REML) (Ch. 27, Lynch & Walsh, 1998), or a prediction error method (Ergon, 2022a(Ergon, , 2022b. Note, however, that REML equations for a single generation will be indeterminate, that is, we can find r 1,t and r 2,t , but not all the residual covariances in the expressions for r 1,t and r 2,t . For use in a general multivariate BLUP equation, we also need the n × n additive genetic relationship matrix A t (Ch. 26, Lynch & Walsh, 1998), and the For comparisons with predictions based on selection gradients (GRAD), as in Equation (5), we also need the phenotypic covariance , with diag([ • ]) denoting diagonal matrices. Also P is assumed to be constant.
From the mixed model according to Equation (7) follows the multivariate BLUP equation in matrix form (Henderson, 1950;Ch. 26, Lynch & Walsh, 1998), Note that the derivation of Equation (8) does not necessarily require an assumption of normal data (Robinson, 1991).

| Updating of mean reaction norm parameter values
When Equation (8) is applied on any given parent generation, the expected mean values of the vector elements in x t will be zero, that However, owing to different fitness (number of descendants) among the individuals in the parent generation, the corresponding mean values in the offspring generation before new reproduction will be different from â′ 1,t , etc., and these within-generation differences may be used for updating of the mean reaction norm parameters a 1,t , a 2,t , b 11,t , b 12,t , b 21,t , and b 22,t in Equations (6a) and (6b). After this updating, the offspring are ready to become new parents, again with

Selection will thus result in within-generation incremental changes
in the mean values of the estimated random effects from parents to offspring before reproduction, generally given by Equation (1), is the corresponding mean value for the descendants of individual i in generation t. In Equation (9), z parents t stands for the mean value of any one of the estimated random effects from Equation (8), that is, â′ 1,t , â′ 2,t , b′ 11,t , b′ 12,t , b′ 21,t , or b′ 22,t , while z offspring t is the corresponding mean value for the offspring after selection but before reproduction. These incremental changes should thus at each generation be used for updating of the mean reaction norm parameter values before the offspring become new parents. For this purpose, we need an additional assumption, which follows from Assumption 3 above, because ẑ′ i,t and ẑ′ descendants i,t are estimates of an additive genetic component of a reaction norm parameter, and thus have no nonadditive components: Assumption 7. There are no expected fitness-weighted changes in estimated random effects from individual parents to their descendants after selection but before reproduction, that is, When Assumption 7 is applied on Equation (9), we obtain the Theorem 1. In a population that is adequately described by Equation (7) value a 1,t , a 2,t , b 11,t , b 12,t , b 21,t , or b 22,t , while ẑ′ i,t is the corresponding estimated individual value â′ 1,i,t , â′ 2,i,t , b′ 11,i,t , b′ 12,i,t , b′ 21,i,t , or b′ 22,i,t , in the random effects vector x t in Equation (8).
From Theorem 1 follows the incremental changes in mean traits according to

With correct initial mean reaction norm parameter values, this
gives y 1,t and y 2,t according to Equation (6b), and this can be generalized to higher dimensions, with p > 2 and q > 2.
For finite population sizes, there will be drift in the mean reaction norm parameter values, owing to random errors in the covariance computations according to Equation (11). However, as we will see in the simulations, there are also other sources of drift.
From the BLUP theory above follows the following theorem: Theorem 2. For a reaction norm evolutionary system with p = 2 phenotypic traits and q = 2 environmental cues, Equation (8)  See Appendix A for proof, and results in Section 3.
For comparisons with the multivariate breeder's equation, the results from Equation (12) can also be found by a generalization of the GRAD incremental parameter changes according to Equation (5):

| 7 of 15
ERGON Theorem 3. Equation (12) can be reformulated as an extension of Equation (5), See Appendix B for proof, and simulation results in Section 3, where Equations (12) and (13) give identical results for population sizes n ≥ 2. The simulations also show that the results from Equation (13) are close to the results from a corresponding version of the multivariate breeder's equation (4), with declining differences for increasing population size.

| Example case without plasticity
When all plasticity slope parameter values are zero, Equation (6a) gives y j,t = a j,t , and the individual phenotypic traits For p traits, this leads to Equation (8) with A proof of Theorem 4 is given in Ergon (2022c), although then relying on a comparison with the multivariate breeder's equation.
Here, Theorem 1 makes it into an independent proof. Note that in this case Equation (13) degenerates into Equation (2).

| Errors in estimated random effects variances
The fact that updated mean reaction norm parameter values, and thus also updated mean phenotypic trait values, are found by the use of Robertson's secondary theorem of natural selection applied on estimated random effects, as stated in Theorem 1, quite generally shows that the variances of these effects are underestimated. This is easily seen in the case without plasticity, given by Equation (14), in which case Equation (2) can be formulated as which should be compared with the result following from Theorem 1, The underestimation of the variances G jj = E a �2 j,i,t becomes especially transparent with diagonal G and P matrices, where we from Equations (15a) and (15b), with the use of 2 In this case, we thus find that is, var � a � j,i,t = � G jj < G jj for 2 j > 0. See Ergon (2022c) for simulation results.

| Adjustments for overlapping generations
With surviving parents, only a fraction f t < 1 of a given generation are offspring from the previous generation. The incremental changes in mean reaction norm parameter values from one generation to the next are then reduced accordingly, and by the use of Equation (11), we thus obtain where z t is any one of the mean reaction norm parameters, while ẑ′ i,t is the corresponding estimated individual random effect. As verified in the simulations, this will slow down responses on environmental changes, with reduced mean fitness as consequence. (13) 3 | S IMUL ATIONS

| The aim of the simulations
The aim of the simulations is to verify the theoretical BLUP results by means of a toy example, and the purpose is fourfold. First, it is verified that mean reaction norm parameter values can be updated from generation to generation by means of Robertson's secondary theorem of natural selection (Theorem 1). Second, it is shown that it is possible to disentangle the microevolutionary and plasticity components of for example climate change acclimations as shown in Equation (11) in general, and in Equation (12) for the special case with A t = I n . Third, it is verified that the dynamical BLUP and GRAD results for the incremental changes in mean reaction norm parameter values are identical for population sizes n ≥ 2, provided that A t = I n (Theorems 2 and 3).
Fourth, it is shown that the GRAD results are erroneous for populations with genetic relatedness between the individuals, that is, for A t ≠ I n .

| Description of toy example
In a toy example in Ergon (2022a) Here, the example is extended to include a second environmental variable with a noisy positive trend in mean value, and with varia- tions from year to year that are somewhat positively correlated with the variations in spring temperature. This input may for example be a measure of spring rainfall. The example also includes a second adaptive phenotype, which might be the breeding habitat, as discussed in Chalfoun and Schmidt (2012). We thus have a microevolutionary system with two environmental cues and two phenotypic traits, similar to the theoretical example case.
The individual (mid-parent) fitness values are integers from 0 to 4, with number of descendants as unit, and cases with both nonoverlapping and overlapping generations are simulated. The population size is assumed to be constant, which implies that not all descendants survive until reproduction. A constant population size is not essential for the principal results, but it simplifies the simulations.
In the simulations, the two environmental reference values are assumed to be known from historical data, that is, it is assumed that the population was fully adapted to the stationary stochastic environment before the onset of anthropogenic and global climate change around 1970.

| Environmental inputs and reaction norm model
Assume a population that is fully adapted to a stationary stochastic en- where the mean values U 1 ,t and U 2 ,t are ramp functions as shown in Figure 1, while u 1,s,t and u 2,s,t are zero mean and white random variables, that is, without autocorrelation. In a corresponding way assume that 1,t = Θ 1 ,t + 1,s,t and 2,t = Θ 2 ,t + 2,s,t , where Θ 1 ,t and Θ 2 ,t are ramp functions as shown in Figure 1, while 1,s,t and 2,s,t are zero mean and white random variables.

| Fitness function and initial mean reaction norm values
The individual fitness function is assumed to be rounded values of where 1,t and 2,t are the phenotypic values that maximize fitness, while 2 = 10. The discrete values of W i,t (number of descendants) are thus integers from 0 to 4.
In the simulations it is essential that the mean reaction norm parameters are given correct initial values at generation t = 1. We will assume that the phenotypic values are scaled such that the initial mean intercept values in a stationary stochastic environment are a 1,t = a 2,t = 0, and that the initial mean reaction norm slope values are the optimal values in a stationary stochastic environment. These optimal values are the ones that maximize the expected individual fitness according to Equation (18), in a stationary stochastic environment, and thus minimize the criterion functions J 1 = E y 1,i,t − 1,t 2 and J 2 = E y 2,i,t − 2,t 2 . With E y 1,i,t = 0, and substituting E y 2 1,i,t = b 2 11,t 2 U 1,s and E y 1,i,t 1,t = b 11,t U 1,s Θ 1,s , we find we thus find the optimal mean slope value b 11,opt = U 1,s Θ 1,s ∕ 2 In the same way, we find that b 22,opt = −0.5 will minimize J 2 .

F I G U R E 2
Simulation results with population size n = 100, and an additive genetic relationship matrix A t = I n . Mean trait values y 1,t and y 2,t are shown by solid blue lines. Mean reaction norm parameter values are shown by solid green lines (BLUP), dashed blue lines (GRAD), and dotted magenta lines (the multivariate breeder's equation).

| Simulation results
Simulation results with population size n = 100, and additive genetic relationship matrix A t = I n , are shown in Figure 2. New additive genetic (random) effects, and nonadditive effects (residuals), were at each generation drawn from normal distributions in accordance with the given G and P matrices. The BLUP results given by Equation (12) (green lines), and the GRAD results given by Equation (13) (dashed blue lines), are identical for population sizes n ≥ 2. Results given by the multivariate breeder's equation (4) (dotted magenta lines), are somewhat different from the BLUP results. These differences are clearly smaller for a population size of n = 1000. Note that y 1,t and y 2,t lag behind 1,t and 2,t , as shown in Figure 1, which is typical for ramp responses from dynamical systems with time constants.
As a test, the simulations were repeated with a constant additive genetic relationship matrix for a population with a high degree of relatedness among individuals, as shown for n = 6 in Equation (19), Figure 3, this gave different results for the BLUP and GRAD methods. In this case new additive genetic (random) effects at each new generation were found as z t = M √ Az 0,t , where z t stands for a ′ 1,t , a ′ 2,t , b ′ 11,t , b ′ 12,t , b ′ 21,t , or b ′ 22,t , and where the different data vectors z 0,t were drawn from normal distributions in accordance with the given G matrix. Here,

M √
A is the matrix square root, that is, This reduced the variances of the random effects to approximately 75 % of the nominal values. The random effects were also mean centred.
The simulations with results as in Figure 3 were repeated for a case with surviving parents, where only a fraction f t = 0.5 of any given generation are offspring from the previous generation, that is, with the use of Equation (16). As seen in Figure 4, panels (a) and (c), and as must be expected, the responses are slowed down. As can also be seen in panels (a) and (c), the fraction f t < 1 causes the mean trait values y 1,t and y 2,t to lag further behind the phenotypic values 1,t and 2,t that maximize fitness. The result of this is lower mean fitness values, as can be seen in the (identical) plots in panels (b) and (d).
The changes in mean values over 60 generations in Figures 2 and 3, will to some degree be caused by drift. As a test, the total change in a 1,t over 60 generations in a stationary stochastic environment as before t = 10 in Figure 1, was computed. For population sizes from n = 10 to 400, and based on 100 repeated simulations this change had a mean value of approximately zero, and a standard error of approximately 0.1, that is, around 10% of the corresponding changes in panel (a) in Figures 2 and 3. With n = 2, this standard error due to drift increased to 15%. There were

F I G U R E 3
Simulation results with population size n = 100, and an additive genetic relationship matrix A according to Equation ( no noticeable differences between cases with A t ≠ I n and A t = I n .
In order to find an explanation of the drift, the environmental cue u 1,t over 60 generations in a stationary environment was approxi-

| SUMMARY AND D ISCUSS I ON
As shown in Section 2, a general reaction norm model with p phenotypic traits and q environmental cues, that is, with p(1 + q) mean reaction norm parameters, can be formulated as a linear mixed model with fixed and random effects. In this model, the fixed effects will be the mean trait values, while the random effects will be the addi-  (7) could be ordered differently, for example as.
, but then also the G, P, U t and Z t matrices must be reorganized accordingly.
The development of the dynamical BLUP model in Equation (8) relies on Assumptions 1, 2, 3, and 7, as given in Section 2, while Although the development of the dynamical BLUP model in Equation (8) does not rely on normal data (Robinson, 1991), it still follows from Theorems 2 and 3 that the BLUP results will be incorrect in cases with non-normal data and an additive genetic relationship matrix A t = I n . A natural conclusion is therefore that non-normal data will give erroneous BLUP results also when A t ≠ I n , at least for the dynamical BLUP model used here.
In the derivation of the linear mixed model in Equation (7), I assumed linear reaction norms as functions of environmental cues u 1,t , u 2,t , etc., but there is nothing in the theory that prevents us from the use of nonlinear reaction norms that are also functions of u 1,t , u 2,t , u 2 1,t , u 2 2,t , etc., as discussed in Gavrilets and Scheiner (1993) and Ergon (2018).

F I G U R E 4
Simulation results with population size n = 100, and additive genetic relationship matrix A t ≠ I n according to Equation (19). In panels (a) and (c), mean trait values y 1,t and y 2,t and mean reaction norm parameter values a 1,t and a 2,t for f t = 1 are shown by solid blue and solid green lines, respectively. Dashed blue and magenta lines show the corresponding responses with a fraction f t = 0.5 of new offspring in the population at all generations. Phenotypic values 1,t and 2,t that maximize fitness are added as weak dotted blue lines. Panels (b) and (d) show identical plots of mean fitness for f t = 1 (solid green lines) and f t = 0.5 (dashed magenta lines).
As shown by Theorem 1, updating of the mean reaction norm parameter values from generation to generation can be done by applying Robertson's secondary theorem of natural selection on the vector elements â′ T 1,t , etc., in the estimated random effects. This shows that the well-known underestimation of the variances of the random effects, is just what is needed, in order to find the correct incremental changes in mean reaction norm parameter values. The resulting dynamics will depend on the additive genetic and residual covariance matrices G and R t , as well as of the additive genetic relationship matrix A t for the parent generation. It is worth noticing that the additive genetic relationship matrix for the offspring generation does not affect the updating of mean reaction norm parameter values. This is shown in the theoretical derivations, but it is also a natural consequence of the fact that the dynamical BLUP model, just as the multivariate breeder's equation, is an evolutionary state-space model (Ergon, 2018).
For cases with A t ≠ I n , it should be noted that A t is included in the BLUP matrix equation (8) via G t , just as in the standard equation used in domestic breeding (Ch. 26, Lynch & Walsh, 1998). It is also reassuring to know that the BLUP equation asymptotically, for A t → I n , gives results that are identical to the GRAD results based on the multivariate breeder's equation (Theorems 2 and 3). Note that results from the use of the multivariate breeder's equation are equal to the GRAD results only asymptotically, that is, for population size n → ∞. For cases without plasticity, and with A t = I n , the results from all three methods are identical, independent of population size (Theorem 4).
Theorems 1, 2, and 3 were verified in simulations with the use of A t = I n , and population sizes down to n = 2, while Theorem 4 was verified by simulations in Ergon (2022c). Simulation results with A t ≠ I n were found by the use of a constant and possibly unrealistic additive genetic relationship matrix, but these results still serve the purpose of showing that the BLUP and GRAD results with A t ≠ I n are different.
The changes in mean trait and mean reaction norm parameter values over 60 generations that can be seen in Figures 2 and 3, are partly caused by drift. As indicated at the end of Section 3, the major sources of drift are the stochastic nature of the input environmental variables u 1,t and u 2,t , as well as of the phenotypic values 1,t and 2,t that maximize fitness. In the simulations, the population size plays a vital role only when n < 10.
In the simulations, the G and P parameters are assumed to be known, although they are normally not, and this is the case also for the reference environment and the initial mean reaction norm parameters. However, as shown in Ergon (2022aErgon ( , 2022b, these parameters can be found by a prediction error method (PEM) using all available input-output experimental or field data, including individual fitness data. Ergon (2022aErgon ( , 2022b combined PEM with the GRAD model, that is, with extensions of Equation (5), but as will be reported separately, PEM works just as well when combined with the BLUP model in Equation (8). BLUP/PEM will in fact be much better than GRAD/ PEM when it comes to the identification of the reference environment. As seen for GRAD/PEM in Ergon (2022aErgon ( , 2022b, and as will be reported separately for BLUP/PEM, the essential feature of PEM in the present setting is to predict the reaction norms as well as possible, while the additive genetic and environmental parameter values may be very much influenced by random measurement errors and modeling errors.
For the interested reader, the essential steps in the proposed dynamical BLUP method are summarized in Appendix C, where also the procedure for PEM system identification from laboratory or field data is included.
Parameter estimation by means of restricted maximum likelihood (REML) applied on data from a single generation is not possible for plastic organisms. A simple reason for this is that the variance components r j,t in the diagonal residual covariance matrix R t are functions of several unknown variances of the nonadditive effects in the reaction norm model, that is, r 1,t , r 2,t , etc. can be found, but not all the residual covariances involved. The reason behind this problem is that there are confounding variables in the reaction norm equations. In for example Equation (3), all the terms in a � i,t + b � i,t u t + i,t + i,t u t are confounded, which implies that the use of REML at a single generation will not give details of the G and R t matrices. As pointed out above, a solution is to use a prediction error method (PEM) where data from all generations are used. Note that this confounding is not a problem in the modeling process, where a t and b t in a � i,t = a i,t − a t and b � i,t = b i,t − b t, respectively, are found from Equation (11), and where a ′ i,t , b ′ i,t , i,t and i,t are samples with given variances.
With surviving parents, only a fraction f t of the population will be offspring from the previous generation, and as shown in Equation (16), the incremental changes in mean reaction norm parameter values will then be reduced accordingly. As verified in simulations, and as must be expected, the responses on environmental changes will then be slowed down, with reduced mean fitness as consequence. Note that the last term on the righthand side of Equation (1), the Price equation, in general, is affected by surviving parents, but that this term according to Assumption 7 does not affect the dynamical BLUP updating from generation to generation.

ACK N OWLED G M ENTS
The GRAD model in Equation (5) was inspired by a comment by Michael B. Morrissay in a review of a quite different manuscript. I thank Bruce Walsh and Allen Moore for constructive review comments, and Tormod Ådnøy for input regarding matrix details in BLUP modeling. I thank the University of South-Eastern Norway for its support and funding.

CO N FLI C T O F I NTE R E S T S TATE M E NT
There are no competing interests.

DATA AVA I L A B I L I T Y S TAT E M E N T
MATLAB code for simulations is given in Supporting information archived on bioRxiv doi: https://doi.org/10.1101/2023.04.09.536146.

A PPE N D I X A
Proof of Theorem 2.
Lemma 1. Assuming an additive genetic relationship matrix A t = I n , the mean values of the estimated fixed effects in Equation (8), for j from 1 to p, are ŷ j,t = y j,t , while the mean values of the estimated random effects, for j from 1 to p and k from 1 to q , are â� j,t = 0 and b� jk,t = 0.
Proof. Equation (8)   Having established that A t = I n results in â� j,t = 0, b� jk,t = 0 and ŷ j,t = y j,t , it is time to find how fitness affects the incremental changes in mean reaction norm parameter values. For clarity of presentation, we here limit the number of phenotypic traits to p = 2, and the number of environmental cues to q = 2, as we also did in parts of Section 2.
We also introduce the diagonal n × n matrix F t = diag w i,t , and the following two lemmas: Lemma 2. The mean value of the product of F t and an estimated random effect vector ẑ t in Equation (8), where ẑ t =â � j,t for j from 1 to p, or ẑ t =b � jk,t for j from 1 to p and k from 1 to q, is F tẑt = Δz t , that is, the incremental change in the mean reaction norm parameter value z t .
Lemma 3. The mean value of the sum F t y j,t − F t 1 n y j,t , for j from 1 to p, is cov w i,t , y j,i,t .
Proof. F t y j,t − F t 1 n y j,t = 1 n ∑ n i=1 w i,t y j,i,t − w t y j,t = cov � w i,t , y j,i,t � .
Assuming p = 2 and q = 2, for simplicity, and using that ŷ j,t = y j,t (Lemma 1), multiplication of Equation (A2a) from the left with the block diagonal matrix F t = I 6 ⨂ F t , and insertion of F −1 tFt between Z T t R −1 t Z t + G −1 and ⨂ I nxt , we find by the use of the mixed model Since all elements in Z T t R −1 t Z t + G −1 and Z T t R −1 t are multiplied by F t , Equation (A3a) can be reformulated as Taking mean values of all vector elements, will according to Lemma 2 and 3 finally give (A1a) a � j,t + u 1,tb � j1,t + u 2,tb � j2,t + … + u r,tb � jq,t = y j,t − 1 nŷj,t , (A1b) a � j,t + u 1,tb � j1,t + u 2,tb � j2,t + … + u r,tb � jq,t = y j,t −ŷ j,t , (A2d) Q jxt = y j,t −ŷ j,t .